function[IRF,Phi,e,A0]=VAR_model(data,hor,p,identification,scale)
% Viet Hoang Dinh
% Monash University
% Edited December 2024

Y=data;
[~,rhs,lhs]=lag_variable(Y,p);

[T N]=size(lhs);
M=size(rhs,2);
Phi=(rhs'*rhs)\(rhs'*lhs);
e=lhs-rhs*Phi;
sigma=e'*e/(T-M);

Phi_int=Phi(2:end,:);
F_var = [Phi_int';eye(N*(p-1)) zeros(N*(p-1),N)];
bigeye= [eye(N) zeros(N,N*(p-1))];

%% Identification
% identification=1 if Cholesky
% identification=2 if Long-Run

if identification==1 % Cholesky
    if scale==0
     A0=chol(sigma,'lower');
    elseif scale==1
       A0=chol(sigma,'lower');
       for i=1:N
           A0(:,i)=A0(:,i)./A0(i,i);
       end  
    end
elseif identification==2 % Long-run
   invC =  eye(p*N)-F_var;
      C =  inv(invC);
      C =  C(1:N,1:N);
      D =  chol(C*sigma*C','lower');
     A0 = inv(C)*D;      
end

%% IRF estimation
 % IRFs
for j= 1:hor+1
    IRF(:,j) = reshape(bigeye*(F_var^(j-1))*bigeye'*A0,[],1);
end


end